---
title: "Vignette"
author: "Dai Yamao"
date: "`r format(Sys.time(), '%Y-%m-%d')`"
output:
  html_document: default
  pdf_document: default
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE, fig.width = 10, fig.height = 5)
rm(list=ls())

require(memisc)
require(stringi)
library(stargazer)
library(coefplot)
library(interplot)
library(sjPlot)
library(sjmisc)
library(sjlabelled)
library(stargazer)
library(stringr)
library(tidyverse)
library(ggplot2)
library(psych)
library(margins)
library(dplyr)
library(modelsummary)
library(makedummies)
library(ggeffects)
library(estimatr)
library(ggpubr)
library(BalanceR)

library(haven)
library(summarytools)
library(tidyverse)
library(ggplot2)

dat <- read_csv("data/PSL_Iraq2025.csv")
dat <- dat[-1,]

dat$Q17a[is.na(dat$Q17a)] <- 0
dat$Q17b[is.na(dat$Q17b)] <- 0
dat$Q17c[is.na(dat$Q17c)] <- 0
dat$Q17d[is.na(dat$Q17d)] <- 0

dat <- dat |>
  mutate(Support = Q17a + Q17b + Q17c + Q17d)
dat$Support[dat$Support == 0] <- NA

dat <- dat |>
  mutate(Treatment = case_when(Q17a >= 1 ~ "Scrapped",
                               Q17b >= 1 ~ "Separate_law",
                               Q17c >= 1 ~ "Age_limit",
                               Q17d >= 1 ~ "Obliged"))
dat$Treatment <- as.factor(dat$Treatment) 

dat <- dat |>
  mutate(Gender = D1,
         Age = D2,
         Education = D3,
         Income = D5)

dat <- dat |>
  mutate(sunni = ifelse(D6 == 1, 1, 0),
         shia = ifelse(D6 == 2, 1, 0),
         christian = ifelse(D6 == 3, 1, 0),
         kurd = ifelse(D6 == 4, 1, 0),
         others = ifelse(D6 == 5, 1, 0))

dat <- dat |>
  mutate(Sect = dplyr::recode(D6, 
                       "1" = "Sunni",
                       "2" = "Shia",
                       "3" = "Christian",
                       "4" = "Kurd",
                       "5" = "others"))

# Islamism
dat <- dat |> 
  mutate(Islamism = case_when(Q2_1 >= 67 ~ "High",
                              Q2_1 <= 66 & Q2_1 >= 34 ~ "Middle",
                              Q2_1 <= 33 ~ "Low"))

# Tribalism
dat <- dat |> 
  mutate(Tribalism = case_when(Q2_3 >= 67 ~ "High",
                               Q2_3 <= 66 & Q2_3 >= 34 ~ "Middle",
                               Q2_3 <= 33 ~ "Low"))
dat$Tribalism <- as.factor(dat$Tribalism)

# patronage_religious_daily_necessities
dat <- dat |> 
  mutate(ptr_rel_daily = case_when(Q4 >= 4 ~ "High",
                                   Q4 <= 2 ~ "Low",
                                   Q4 == 3 ~ "other"))
dat$ptr_rel_daily <- as.factor(dat$ptr_rel_daily)

# patronage_Tribe_daily_necessities
dat <- dat |> 
  mutate(ptr_trb_daily = case_when(Q6 >= 4 ~ "High",
                                   Q6 <= 2 ~ "Low",
                                   Q6 == 3 ~ "other"))
dat$ptr_trb_daily <- as.factor(dat$ptr_trb_daily)

# Shia parties supporter
dat <- dat |> 
  mutate(Shia_party_supporter = (Q4_1 + Q4_2 + Q4_3)/3) |>
  mutate(Shia_party_supporter = ifelse(Shia_party_supporter >= 7, "Shia Party Supporter", "others"))
dat$Shia_party_supporter <- as.factor(dat$Shia_party_supporter)
```

# balance check
```{r}
balance <- BalanceR(data = dat, 
                    group = Treatment,
                    cov = c(Gender, Age, Education, Income))
print(balance, digits = 3)
summary(balance, digits = 5)
plot(balance, point.size = 5, text.size = 18)
```

# boxplot
```{r}
dat |> 
  ggplot() +
  geom_boxplot(aes(x = Treatment, y = Support)) +
  scale_x_discrete(na.translate = FALSE)
```

# regression
```{r}
reg01 <- lm(formula = Support ~ Treatment,
            data = dat)
summary(reg01)

reg02 <- lm(formula = Support ~ Treatment +
              Gender + Age + Education + Income,
            data = dat)
summary(reg02)

reg03 <- lm(formula = Support ~ Treatment * Islamism +
              Gender + Age + Education + Income,
            data = dat)
summary(reg03)

reg04 <- lm(formula = Support ~ Treatment * Shia_party_supporter +
              Gender + Age + Education + Income,
            data = dat)
summary(reg04)

reg05 <- lm(formula = Support ~ Treatment * ptr_rel_daily +
              Gender + Age + Education + Income,
            data = dat)
summary(reg05)

reg06 <- lm(formula = Support ~ Treatment * Tribalism +
              Gender + Age + Education + Income,
            data = dat)
summary(reg06)

reg07 <- lm(formula = Support ~ Treatment * ptr_trb_daily +
              Gender + Age + Education + Income,
            data = dat)
summary(reg07)


#
reg13 <- lm(formula = Support ~ Treatment * Q2_1 +
              Gender + Age + Education + Income,
            data = dat)
summary(reg13)
```

# regression table
```{r}
reg_table <- mtable("Model 1" = reg01, "Model 2" = reg02, 
                    "Model 3" = reg03, "Model 4" = reg04,
                    "Model 5" = reg05, "Model 6" = reg06,
                    "Model 7" = reg07,
               summary.stats=c("adj. R-squared","F","p", "N"))
print(reg_table)
write_html(reg_table, file = "result/result_vignett.html")
```


# marginal effect
```{r}
plot_model(reg01, type = "pred", terms = c("Treatment"),
           axis.labels = "Treatment", title = "")
plot_model(reg02, type = "pred", terms = c("Treatment"),
           axis.labels = "Treatment", title = "")

plot_model(reg03, type = "pred", terms = c("Treatment", "Islamism"),
           axis.labels = "Treatment", title = "")
plot_model(reg13, type = "pred", terms = c("Q2_1" ,"Treatment"),
           axis.labels = "Treatment", title = "")

plot_model(reg04, type = "pred", terms = c("Treatment", "Shia_party_supporter"),
           axis.labels = "Treatment", title = "")

plot_model(reg05, type = "pred", terms = c("Treatment", "ptr_rel_daily"),
           axis.labels = "Treatment", title = "")

plot_model(reg06, type = "pred", terms = c("Treatment", "Tribalism"),
           axis.labels = "Treatment", title = "")

plot_model(reg07, type = "pred", terms = c("Treatment", "ptr_trb_daily"),
           axis.labels = "Treatment", title = "")

```

# prots
```{r}
p1 <- plot_model(reg01, type = "pred", terms = c("Treatment"),
           axis.labels = "Treatment", title = "")

p2 <- plot_model(reg03, type = "pred", terms = c("Treatment", "Islamism"),
           axis.labels = "Treatment", title = "")

p3 <- ggarrange(p1, 
          p2 + theme(axis.title.y = element_blank()), 
          nrow = 1, ncol = 2, align = "hv")
p3

ggsave(filename = "result/scenario.png",
       plot = p3,
       width = 9,
       height = 4,
       units = "in", 
       dpi = 600,
       bg = "white")
```
